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ABSTRACT 


The  development  of  a  two-dimensional  computer  code  to  simulate 
detonation  transfer  between  explosive  layers  is  described.  The  code  is 
based  on  previous  models  developed  at  the  Laboratory  for  Computational 
Physics  and  Fluid  Dynamics  at  Naval  Research  Laboratory,  Washington,  DC 
to  study  the  structure  of  layered  detonations  and  the  details  of  detonation 
transmission  from  one  medium  to  another.  The  code  has  been  configured 
to  simulate  experiments  conducted  at  the  University  of  Michigan  on  the 
lateral  transfer  of  detonation  and  shock  phenomena  between  different 
gaseous  layers.  Preliminary  calculations  with  the  code  show  that  the 
computations  produce  many  of  the  structures  seen  in  the  Michigan 
experiments  and  also  provide  detailed  descriptions  of  the  detonation 
transmission  and  evolving  structure. 
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NUMERICAL  SIMULATION  OF  DETONATION  TRANSFER 
BETWEEN  GASEOUS  EXPLOSIVE  LAYERS 


1.  INTRODUCTION 


The  transmission  and  reflection  of  shock  and  detonation  waves  through  layered  materials 
consisting  of  combinations  of  inert  materials  of  different  shock  impedences,  or 
combinations  of  inert  and  detonable  materials,  are  of  fundamental  importance  in  many 
areas  of  military  interest.  Recent  examples  at  MRL  which  consider  shock  transmission 
normal  to  the  layers  have  included  studies  of  shock  propagation  through 
copper/kapton/copper  layers  for  the  analysis  of  a  shock  compression  conduction  switch  [11, 
and  the  transmission  of  impact  shocks  through  metal /explosive  layers  for  the  analysis  of  the 
jet  initiation  of  covered  explosives  [21.  When  layered  explosives  are  involved  an  interesting 
question  is  whether  a  suitable  combination  of  explosives  can  enhance  the  effectiveness  of 
the  energy  transmission  compared  to  a  pure  explosive  with  the  same  total  energy. 

Interesting  effects  also  occur  when  a  detonation  propagates  along  an  explosive 
layer  and  then  comes  into  contact  with  a  bounding  explosive.  If  the  inner  and  bounding 
explosives  are  identical,  the  layered  detonation  pattern  represents  the  diffraction  that 
occurs  when  a  detonation  propagates  past  a  step  or  an  increase  in  cross-sectional  area. 

With  a  suitable  choice  of  symmetry  axis  this  system  can  also  represent  the  propagation  of  a 
detonation  from  a  smaller  into  a  larger  tube,  or  an  inner  explosive  sheathed  in  an  outer 
explosive. 


When  solid  explosives  are  used  to  study  detonation  transfer  between  layers  it  is 
very  difficult  to  observe  the  complex  interactions  which  develop.  The  problem  is  very 
much  simpler  for  gaseous  detonation  transfer  however,  as  pulsed  laser  Sehlieren 
photography  is  an  established  method  for  studying  complex  structures  in  gaseous 
explosives.  At  the  University  of  Michigan  this  technique  has  been  combined  with  a 
specially  designed  shock  tube  to  visualize  transmitted  shock  and  detonation  structures  in 
layered  detonations  for  an  extensive  combination  of  explosives  [3-61.  The  experimental 
arrangement  consists  of  two  adjacent  1.6  cm  x  1.6  cm  square  detonation  tubes  over  three 
metres  in  length.  A  test  section  of  length  16.2  cm  or  19  cm  is  included  and  the  different 
gases  are  initially  separated  by  a  nitrocellulose  film  of  about  60  n.n  thickness.  A  stable 
detonation  is  initiated  in  the  primary  detonation  tube  and  the  interaction  between  the  two 
mixtures  begins  when  the  detonation  in  the  primary  tube  comes  into  contact  with  the 
film.  In  the  initial  phase  of  the  interaction  an  explosive  bubble  or  blast  wave  propagates 
into  the  secondary  mixture.  Further  progress  of  the  interaction  then  depends  on  the 
properties  of  both  of  the  explosive  mixtures. 


In  a  series  of  experiments  recently  reported  [51,  the  primary  detonation  tube 
contained  a  stoichiometric  mixture  of  H2  and  O2  while  the  secondary  tube  contained  a 
variable  mixture  of  Ho  and  02-  The  exact  composition  of  the  explosive  in  the  secondary 
tube  is  described  by  the  equivalence  ratio  $ ,  which  is  defined  to  be  the  ratio  of  the  amount 
of  fuel  to  oxygen  for  the  given  system,  divided  by  the  same  ratio  for  a  stoichiometric 
system.  In  the  experiments  described  in  131  the  equivalence  ratio  of  the  mixture  varied 
between  0.15  and  3.5.  Three  main  modes  of  interaction  were  observed  as  the  equivalence 
ratio  of  the  bounding  explosive  was  increased.  In  the  first  mode,  (<p  <  0.25),  only  an 
oblique  shock  was  induced  in  the  mixture  and  there  was  no  detonation.  In  the  second  mode 
(0.35  <  <p  <  1.0)  detonation  was  initiated  directly,  while  in  the  third  mode 
(0.35  <  1 p  <  4.5)  detonation  was  initiated  only  after  a  Mach  reflection  of  the  induced 
oblique  shock  from  the  lower  wall. 

Shock  polar  analysis  of  the  steady  state  interaction  and  the  reflection  from  the 
lower  wall  have  been  used  to  analyse  the  different  modes  of  interaction.  A  two-gamma 
method  described  by  Liou  [71  has  been  used  to  compute  the  oblique  shock  and  detonation 
angles  in  the  bounding  explosive.  The  detonation  angles  computed  using  this  method  follow 
the  trend  of  the  experimental  results,  although  the  computed  results  are  significantly  below 
the  measured  values.  Accurate  analysis  of  the  transient  processes  which  occur  as  the 
explosive  bubble  first  propagates  into  the  bounding  explosive  cannot  be  treated  using  this 
approach  however  and  require  detailed  numerical  simulation. 

The  purpose  of  this  report  is  to  describe  the  development  of  a  two-species 
reactive  hydrocode  to  model  these  transient  processes.  The  code  described  here  is  based 
on  the  Reactive  Shock  Models  described  by  Oran  and  Boris  [81.  The  convective  transport 
equations  are  solved  using  the  Flux-Corrected  Transport  (FCT)  algorithm  developed  in  the 
Laboratory  for  Computational  Physics  and  Fluid  Dynamics  at  the  Naval  Research 
Laboratr ;  y  19,10,111.  The  particular  FCT  algorithm  used  has  fourth-order  accuracy  and 
can  reproduce  discontinuities  with  small  dispersion  errors  and  minimal  numerical 
diffusion.  It  has  already  been  used  successfully  to  study  Mach  reflection  of  detonation 
waves  in  a  variety  of  simulations  (12,131. 

As  the  basic  reactive  shock  model  used  in  this  work  has  been  fairly  extensively 
described  elsewhere  [12,141  our  description  of  the  two-species  model  (Section  2),  and  its 
numerical  implementation  (Section  3),  is  relatively  brief  and  concentrates  mainly  on  those 
aspects  unique  to  the  geometry  and  two-species  nature  of  the  problem.  A  detailed 
discussion  of  the  testing  at  each  stage  of  code  development  and  a  comparison  of  the 
numerical  simulations  with  the  experimental  results  is  presented  in  Section  4. 


2.  DESCRIPTION  OF  MODEL 


The  code  developed  here  is  based  on  similar  models  described  recently  by  Guirguis  et  al. 
[121.  The  Euler  equations  for  compressible  flow  are  solved  assuming  perfectly  rigid  and 
smooth  walls.  Two-dimensional  rectangular  geometry  is  employed  with  uniform  grid 
spacing  and  fixed  cell  sizes  in  the  x  and  y  directions.  The  complexity  and  extent  of  the 
interactions  behind  the  leading  shock  front  make  the  use  of  moving  or  adaptive  gridding 
techniques  [81  inappropriate  for  this  particular  problem. 

The  equations  for  the  conservation  of  mass,  momentum  and  energy  are 
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Here  p  is  the  density,  u  and  v  the  x  and  y  components  of  the  velocity  vector  V,  P  is  the 
scalar  pressure,  and  E  is  the  total  energy  per  unit  volume,  defined  by 

E  =  pe  +  ^  p(u2  +  v2) ,  (4) 

where  e  is  the  specific  internal  energy. 

Equations  (1)  through  (3)  are  solved  using  the  JPBFCT  algorithm  (a  version  of 
ETBFCT,  documented  in  1151)  and  operator  splitting  in  the  x  and  y  directions.  Hie  rigid 
barrier  separating  the  two  detonation  tubes  runs  along  the  x  direction  and  is  included  in  the 
calculation  during  the  y  integration.  Depending  on  whether  the  x-coordinate  is  before  or 
after  the  end  of  the  barrier,  the  y  integration  is  split  into  either  two  or  one  loops,  with  the 
end  points  of  the  loops  being  the  upper  and  lower  boundaries  of  each  tube  separately,  or  the 
upper  and  lower  boundaries  of  the  combined  tube  (see  Figure  1). 

The  total  density  p  is  the  sum  of  the  densities  of  species  one  and  two 


p  ~  pSPl  +  pSP2'  (5) 

Each  species  is  convected  separately  with  the  fluid  velocity  and  normalised  after  each  time 
step.  The  species  are  initially  separated  in  the  upper  and  lower  detonation  tubes  but  are 
allowed  to  mix  freely  as  the  detonation  in  the  upper  tube  encounters  the  end  of  the  dividing 
barrier. 

A  perfect-gas  equation  of  state  is  used  for  each  species,  ie.  P.  =  p.R.T,  where 
Rj  is  the  specific  gas  constant  for  species  i,  and  T  is  the  equilibrium  temperatuVe'of  both 
species.  In  all  of  the  calculations  described  below,  the  full  details  of  the  chemical 
reactions  are  not  included  in  the  model.  Instead  we  use  an  induction  parameter  model  that 
reproduces  the  essential  features  of  the  chemical  reaction  and  energy  release  process.  In 
this  model,  the  chemical  induction  time  and  energy  release  rate  are  tabulated  as  a  function 
of  temperature,  pressure,  and  stoichiometry.  These  have  been  obtained  by  integrating  the 
full  set  of  elementary  chemical  reactions  for  the  hydrogen-oxygen-argon  system  using  the 
CHEMEQ  integration  code  [161.  Then  a  quantity  called  the  induction  parameter  is  defined 
and  convected  with  the  fluid  in  a  Lagrangian  manner.  This  parameter  records  the 
temperature  history  of  a  fluid  element  and,  when  the  element  has  been  heated  for  long 
enough,  energy  release  is  initiated.  This  model  for  including  the  properties  of  a  chemical 
reaction  mechanism  in  a  numerical  simulation  was  described  originally  by  Oran  et  al.  [201, 
and  has  been  developed  further  by  Kailasanath  et  al.  [211  and  Guirguis  et  al.  [121. 

The  induction  time  t°  ,  defined  to  be  the  time  during  which  no  energy  is 
released,  is  fitted  to  an  expression  of  the  form 


(6) 


i°(T,P)  =  At(Po/P)exp(Er/Kn 


In  the  second  step,  reactants  are  converted  to  products  according  to  the  finite  reaction  rate 


du/dt  =  -  <dArexp(-Er/R  T) ,  (7) 

where  w  is  the  explosive  mass  fraction  for  species  one  or  two.  The  constants  A  ,  E^ ,  Ar 
and  Er  are  calculated  from  the  results  of  integrating  the  full  set  of  chemical  reactions  for 
the  system  [161.  If  f  denotes  the  fraction  of  induction  time  elapsed  at  time  t,  then 

£  -  -2—  ,  (8) 
dt  t°(T,P) 

where  f(0)  =  0.  However,  as  explained  ir  [101,  equation  (8)  is  rewritten  as 
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Hie  final  temperature  and  pressure  within  a  cell  are  given  by 
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where  n .  is  the  fraction  of  species  i  within  the  cell,  and  aEj  ,  Cvj  and  mwj  are  the  energy 
release,  specific  heat  and  molecular  weight  for  species  i. 


3.  NUMERICAL  SOLUTION 


In  a  moving  system,  the  time  derivatives  appearing  in  equations  (6)  and  (9)  denote  a 
substantial  derivative  following  the  flow  of  the  system.  When  combined  with  the 
continuity  equation,  they  can  be  rewritten  in  the  same  form  as  equations  (1)  to 
(3),  i.e. 
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for  each  species  i.  The  solution  procedure  is  then  as  follows:  the  variables 

p,  pu,  pv,  E,  pr). ,  pa.  and  poi.f.  are  advanced  over  the  time  step  fit  using  the  JPBFCT 

algorithm  and  operator  splitting.  1lhe  specific  internal  energy  is  then  calculated, 

the  n-  renormalised,  and  the  explosive  mass  fraction  o>.  is  then  limited  to  a  value  less  than 

or  equal  to  one  using 


a,.  =  min  11.0,  (14) 

The  fraction  u^f .  is  also  limited  according  to 

oi. f.  =  min  (oi.n.,  oi. f.)  (15) 

it  liii 

The  temperature  and  pressure  are  then  calculated  from  equations  (10)  and  (11). 

Next,  a  burn  subroutine  is  called  to  calculate  the  energy  release  if  the  induction 
time  has  elapsed.  First  At,  the  time  remaining  before  the  end  of  the  induction  period,  is 
calculated  and  then  compared  with  the  hydrodynamic  time  step  fit .  There  are  two 
possibilities.  At  is  either  greater  or  less  than  fit .  In  the  former  case  the  fraction  of 
induction  time  elapsed  is  advanced  using  an  explicit  solution  to  the  equation 

3  (  poi .  f .  ) 

-sr-5-  -  PiV'!  (16> 

and  control  returns  to  the  main  program.  In  the  latter  case  the  reaction  variable  oi.  is 
updated  using  an  implicit  solution  to  the  equation  1 
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and  a  new  temperature  and  pressure  are  calculated  using  equations  (10)  and  (11). 

The  hydrodynamic  time  step  fit  is  calculated  using  the  Courant  condition, 


=  CENT  *  min  ( 


I  ve  1 1  +  c 


where  Ai  is  the  cell  length  in  either  the  x  or  y  direction,  vel  is  the  x  or  y  component  of 
velocity  in  each  cell,  and  c  is  the  sound  speed.  In  a  cell  with  a  mixture  of  species  one  and 
two,  the  sound  speed  is  calculated  for  each  species  separately  and  the  greater  of  the  two  is 
used.  In  all  the  calculations  to  be  described  here,  CRNT  has  the  value  0.26. 

A  typical  computational  grid  has  500  cells  in  the  x  direction  and  80  cells  in  the  y 
direction,  with  detonation  initiated  at  the  left  hand  boundary  of  thr  grid  (x  =  0).  The  value 
of  the  detonation  velocity  is  such  that  it  will  typically  take  3,000  to  4,000  cycles  for  the 
detonation  to  reach  the  right  boundary.  With  a  fixed  right  boundary,  this  means  that  in  the 
initial  stages  of  the  calculation,  a  great  deal  of  computational  time  is  wasted  computing 
properties  of  large  areas  of  the  grid  in  which  nothing  is  happening.  This  is  very  time 
consuming  and  hence  expensive.  To  solve  this  problem  we  decided  to  make  the  right 
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boundary  movable.  The  number  of  cells  in  the  z  direction,  NCX,  is  initially  set  to  some 
conveniently  small  value,  say  NCX  =  50,  and  then  a  check  is  made  at  the  end  of  each  time 
step  to  locate  the  right-most  position  of  the  leading  shock.  NCX  is  then  adjusted  to  make 
sure  that  it  is  always  five  cells  ahead  of  the  position  of  the  shock. 

In  each  of  the  calculations  described  in  the  next  section,  the  detonation  in  the 
primary  tube  was  initiated  by  depositing  an  excess  amount  of  energy  into  the  first  ten 
columns  of  the  grid,  thereby  creating  a  blast  wave  that  evolved  into  a  detonation.  The 
amount  of  energy  per  cell  required  to  initiate  detonation  in  the  2D  code  once  the  blast  had 
progressed  beyond  the  end  of  the  barrier  was  typically  an  order  of  magnitude  greater  than 
that  needed  to  initiate  a  stable  CJ  detonation  in  a  one-dimensional  code.  This  means  that 
the  detonation  was  considerably  overdriven  as  it  reached  the  end  of  the  dividing  barrier,  a 
situation  which  did  not  apply  in  the  actual  experiments.  Nevertheless,  once  the  detonation 
passed  the  end  of  the  barrier,  its  velocity  rapidly  approached  the  CJ  value  and  some 
distance  from  the  divider  the  pressure  contours  show  patterns  similar  to  those  observed 
experimentally.  Future  use  of  the  code  however  will  include  an  option  to  input  steady 
profiles  for  a  stable  CJ  detonation  taken  from  a  one-dimensional  code. 


4.  CODE  TESTING  AND  RESULTS 


Most  of  the  experimental  work  at  the  University  of  Michigan  reported  to  date  has  used 
mixtures  of  H2  and  02  in  both  the  primary  and  secondary  detonation  tubes.  As  we  do  not 
yet  have  analytical  expressions  for  induction  times  and  energy  release  times  for  H2  -  02 
mixtures  with  given  equivalence  ratio  <P,  the  results  described  here  use  parameters  for 
stoichiometric  H2  -  O,  diluted  with  Ar  (H2:02:Ar  in  ratio  2:1:7)  which  were  available  from 
previous  simulations  of  shock  tube  experiments  1171.  Values  for  each  of  the  parameters 
defined  in  sections  two  and  three  for  this  mixture  can  be  found  in  Table  1.  A  one- 
dimensional  version  of  the  code  was  used  to  calculate  the  CJ  detonation  velocity  by 
depositing  a  small  amount  of  excess  energy  into  the  system  and  allowing  the  simulation  to 
run  until  the  detonation  velocity  had  reached  a  steady  state.  The  value  found  using  this 
method  was  1600  m/sec. 

The  two-dimensional  code  was  developed  in  several  stages  and  tested  at  the 
completion  of  each  stage.  We  began  with  a  two-dimensional  non-reactive  hydrocode  in 
rectangular  geometry.  The  capabilities  of  the  code  at  this  stage  are  illustrated  in  Figure  2, 
which  shows  pressure  contours  every  400  cycles  after  an  excess  amount  of  energy  was 
deposited  into  a  few  cells  centered  near  the  left  boundary.  Each  of  the  surfaces  is  a  rigid 
reflecting  barrier.  These  figures  clearly  show  the  ability  of  the  operator-split  FCT  code  to 
model  complex  shock  structures.  Mach  reflection  of  the  blast  wave  from  the  top  and 
bottom  boundaries  is  clearly  visible,  and  the  convergent  motion  of  the  triple  points  can  also 
be  seen.  The  interaction  of  the  shocks  reflected  from  the  left  upper  and  lower  comers 
with  the  converging  triple  points  has  produced  quite  a  complicated  shock  structure  after 
1400  cycles.  This  structure  progressively  collapses  onto  itself  and  eventually  results  in  the 
formation  of  a  decaying  plane  wave  moving  to  the  right.  The  contour  at  2600  cycles  shows 
the  shock  structure  just  before  this  stage  is  reached. 

Our  first  modification  was  to  simulate  a  rigid  barrier  dividing  the  left  half  of 
the  grid  into  separate  upper  and  lower  detonation  tubes.  Figure  3  shows  a  sequence  of 
pressure  contours  in  a  calculation  in  which  a  blast  wave  was  initiated  in  the  upper  tube  and 
allowed  to  progress  past  the  end  of  the  barrier  (v.hich  extends  along  half  the  length  of  the 
grid).  The  contour  plot  at  500  cycles  shows  the  blast  wave  completely  confined  within  the 
upper  tube.  The  complex  shock  structure  at  this  stage  is  very  similar  to  the  structure 
shown  in  Figure  2  after  1000  cycles.  After  2600  cycles  the  structure  has  decayed  to  a 
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planar  shock  which  has  progressed  past  the  end  of  the  barrier  and  initiated  a  blast  wave  in 
the  second  tube.  By  3600  cycles  the  shock  is  just  about  to  reflect  from  the  bottom  of  the 
second  tube.  The  progression  of  the  shock  backwards  into  the  second  tube  can  also  be  seen. 

The  next  stage  was  to  make  the  code  reactive  by  including  a  model  that  allows  a 
single  species  to  react,  as  decribed  in  the  previous  sections.  At  this  stage  we  also  set  the 
number  ofcells  in  the  x  and  y  directions  to  600  and  80  respectively,  with  Ax  =  Ay 
=  4.0  x  10-4  m.  TTie  value  of  the  cell  size  was  decided  from  previous  work  on  this 
system.  Figure  4  shows  the  result  of  initiating  a  detonation  in  the  upper  chamber  by 
depositing  an  excess  amount  of  energy  into  the  first  10  columns  of  the  grid.  As  the  code  at 
this  stage  is  a  single  species  code,  the  simulations  represent  experiments  performed  with 
the  same  reactive  mixture  in  both  the  upper  and  lower  tubes.  For  experiments  with  very 
thin  films  sepe rating  the  two  similar  layers,  detonation  is  initiated  almost  instantaneously 
in  the  lower  tube  and  Mach  reflection  of  the  detonation  always  occurs  at  the  lower 
surface.  When  the  film  is  thicker,  a  blast  wave  is  transmitted  and  then  the  secondary 
detonation  is  initiated  behind  the  Mach  reflection  from  the  lower  wall.  Examination  of  the 
pressure  contours  in  Figure  4  shows  similar  behaviour  to  the  experimental  results  found 
when  very  thin  films  are  used,  detonation  is  initiated  directly  by  the  blast  wave  and  then  a 
Mach  reflection  is  formed  on  the  lower  boundary.  The  reflected  wave  patterns  coincide 
nicely  with  the  experimental  results  shown  in  [61. 

The  last  stage  involved  converting  the  program  into  a  two-species  code.  We  did 
this  by  first  defining  new  arrays  and  data  for  a  second  species  and  convecting  the  variables 
as  outlined  in  the  previous  section.  We  then  ran  a  calculation  in  which  the  second  species 
was  treated  as  an  inert  gas.  Figure  6  shows  pressure  contours  up  to  3000  cycles  for  the 
case  where  species  one  is  H2:02:Ar  in  the  ratio  2:1:7.  The  second  species  is  the  same 
mixture  as  the  first,  but  it  is  constrained  so  that  there  is  no  energy  release.  The  contrast 
with  the  results  shown  in  Figure  4  is  quite  interesting,  the  inability  of  the  second  species  to 
release  energy  has  delayed  the  formation  of  the  Mach  stem  and  resulted  in  the  formation  of 
a  complex  reflected  shock  pattern  which  moves  steadily  through  the  system. 

Finally,  we  further  modified  the  code  to  allow  the  second  species  to  detonate 
and  ran  a  simulation  with  the  same  specieB  in  both  tubes  as  a  test  of  the  program.  Pressure 
contours  up  to  2000  cycles  for  this  run  are  shown  in  Figure  6.  These  results  should  be 
compared  with  Figure  4,  which  used  the  same  mixture  in  the  upper  and  lower  chambers  but 
were  generated  using  the  single-species  code.  The  results  are  very  similar  to  one  another 
but  not  quite  identical.  This  could  be  due  to  slight  differences  in  the  initial  condition 
between  the  two  calculations,  or  to  diffusional  effects  between  the  species  across  the 
interface. 


Having  fully  developed  the  code  and  checked  that  it  was  giving  satisfactory 
results,  we  then  ran  a  simulation  with  different  reactive  species  in  the  upper  and  lower 
chambers.  Rather  than  spending  a  considerable  amount  of  time  calculating  induction  times 
and  energy  release  times  for  the  pure  H2-02  mixtures  used  in  the  experiments,  we  decided 
to  retain  the  existing  parameters  for  H2:02:Ar  and  to  vary  those  parameters  in  the  same 
direction  as  the  changes  in  the  corresponding  parameters  for  H2-02  as  p  varied  from  0.125 
to  2.6.  We  ran  the  chemical  kinetics  code  CHEMOD  [81  for  selects!  values  of  initial 
pressure  and  temperature  for  Hg-02  mixtures  with  different  values  of  <p.  Over  the  range 
0.125  to  1.0,  we  found  that  the  induction  time  remained  almost  constant,  while  between  1.0 
and  2.6,  it  more  than  doubled.  The  energy  release  time  decreased  only  slightly  over  the 
full  range.  With  stoichiometric  H2-02  in  the  upper  chamber  and  a  variable  mixture  of  H2 
and  Oo  in  the  lower  chamber,  the  experiments  show  a  change  in  interaction  mode  occurring 
around  p  =  0.25.  Consequently  we  decided  to  change  the  parameter  values  for  our  model 
system  to  simulate  a  p  =  0.25  mixture  in  the  lower  chamber.  As  we  had  found  that  there 
were  only  small  changes  in  both  induction  time  and  energy  release  time 
between  p  -  1.0  and  p  =  0.25,  we  decided  to  leave  these,  and  the  value  for  y, 
unchanged.  We  simulated  the  p  =  0.26  mixture  by  setting  AEg  =  0.070  kcal  and 
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mv?2  =  60.0.  In  doing  this  we  hoped  to  find  Interaction  Mode  1  type  behaviour,  ie.,  the 
secondary  mixture  is  not  initiated  and  only  oblique  shocks  reflecting  regularly  from  the 
lower  wall  are  observed.  The  pressure  contours  for  this  calculation  are  shown  in 
Figure  7.  Direct  initiation  of  the  bounding  explosive  by  the  blast  wave  is  again  observed, 
but  the  reflection  of  this  detonation  from  the  lower  boundary  is  quite  different.  The 
detonation  in  the  bounding  explosive  initiated  by  the  blast  wave  appears  to  be  unstable  and 
is  rapidly  quenched.  The  detonation  is  then  re-initiated  behind  the  shock  reflected  from 
the  lower  wall  and  a  complicated  pattern  of  an  oblique  shock  connected  to  an  oblique 
detonation  appears  to  develop.  This  behaviour  has  been  seen  experimentally  118)  and  would 
be  an  interesting  case  for  further  study. 

We  next  attempted  to  find  Mode  1  behaviour  by  making  further  changes  to  the 
parameters  specifying  species  two  so  that  we  simulated  an  even  leaner  mixture,  ie.,  <p  less 
than  0.25,  but  here  we  ran  into  problems.  While  the  values  for  aE_  and  mwj  chosen  for  <p  = 
0.25  lead  to  a  stable  detonation  in  a  ID  code  with  a  CJ  velocity  or  approximately  1000  m/B, 
no  detonation  in  a  ID  code  could  be  found  when  the  parameter  changes  were  made  in  the 
direction  of  decreasing  <p .  It  appears  that  we  have  run  into  the  detonation  limit  for  our 
model  system.  Remembering  that  the  original  system  is  already  quite  dilute,  this 
behaviour  does  not  seem  unreasonable.  At  this  stage  it  was  decided  that  further  numerical 
simulation  of  the  experimental  results  required  expressions  for  the  induction  time  and 
energy  release  time  for  pure  H2-O2  mixtures.  Calculation  of  these  quantities  is  currently 
in  progress,  and  results  will  be  reported  at  a  later  date. 


5.  DISCUSSION 


The  results  described  in  the  previous  section  are  particularly  encouraging.  When  the 
materials  in  both  tubes  are  identical  the  simulations  show  that  a  detonation  is  initiated 
directly  in  the  second  tube  and  a  Mach  reflection  occurs  at  the  lower  boundary.  When  a 
leaner  mixture  is  used  in  the  second  tube  the  simulations  again  show  structures  very  similar 
to  those  seen  in  experiments.  We  cannot  expect  quantitative  agreement  between  the 
simulations  and  experiments  at  this  stage  as  the  experiments  were  performed  for  mixtures 
of  Hg  and  02,  while  the  simulations  were  performed  for  mixtures  of  H2  and  Oa  diluted  with 
Argon.  These  diluted  mixtures  will  have  a  lower  detonation  velocity  as  well  as  a  more 
curved  detonation  front  due  to  a  longer  induction  time.  Quantitative  comparisons  could  be 
made  if  suitable  constants  could  be  determined  for  the  induction  model  and  energy  release 
rate  expression  described  in  Section  2.  These  can  be  found  by  integrating  a  detailed 
chemical  kinetics  model  [16],  and  checked  by  comparison  with  experimental  values  for 
induction  times  in  undiluted  H2-C>2  mixtures  [71.  Work  along  these  lines  is  proceeding,  and 
we  anticipate  that  quantitative  comparisons  between  the  simulations  and  the  experiments 
will  be  available  at  a  future  date. 

Before  further  calculations  are  made  with  the  code  some  minor  changes  could  be 
made  which  would  facilitate  comparison  of  the  output  with  the  experimental  results. 

Future  runs  should  output  contour  plots  of  density,  temperature,  velocity,  and  reaction 
progress  variables  as  well  as  pressure  contours,  as  these  would  aid  the  interpretation  of  the 
complex  shock  and  detonation  structures.  The  experimental  results  are  produced  by 
Schlieren  photography  and  this  technique  is  sensitive  to  density  gradients,  whereas  the 
results  presented  so  far  have  only  shown  pressure  contours.  Density  contours  would  provide 
a  better  comparison  with  experiment,  and  also  enable  important  leatures  such  as  contact 
surfaces  to  be  seen.  The  experimental  apparatus  also  has  pressure  transducers  mounted  on 
the  top  and  bottom  of  the  detonation  tube  so  that  pressure-time  histories  at  selected 
positions  can  be  obtained.  TV.s  information  could  easily  be  output  from  the  code  so  that 
further  comparison  with  experiment  could  be  made. 
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An  important  change  which  should  be  made  involves  the  method  of  initiation. 
Currently  the  code  is  started  by  depositing  excess  energy  into  the  first  10  cells  in  each  row 
of  the  primary  chamber.  This  creates  a  blast  wave  which  almost  instantaneously  becomes 
an  overdriven  detonation  which  travels  down  the  chamber.  As  the  detonation  reaches  the 
end  of  the  dividing  barrier  it  is  still  highly  overdriven,  in  contrast  to  the  experiment,  but 
soon  decays  to  a  stable  detonation  propagating  at  the  CJ  velocity.  A  much  better  method 
of  initiation,  one  which  parallels  the  experimental  method  much  more  closely,  is  to  use 
output  from  a  ID  reactive  hydrocode  as  the  input  to  the  2D  code.  A  ID  reactive  hydrocode 
was  written  during  the  course  of  this  work  to  check  the  CJ  values  for  the  diluted  H2-O2 
mixtures  used  [191.  This  code  is  also  initiated  by  putting  excess  energy  into  the  first  few 
cells,  but  is  then  allowed  to  run  for  a  sufficient  amount  of  time  so  that  a  constant 
detonation  velocity  is  established.  At  this  stage  variable  profiles  could  be  saved  and  used 
as  input  to  the  2D  code,  which  would  ensure  that  a  stable  CJ  detonation  was  approaching 
the  end  of  the  dividing  barrier  in  the  2D  simulation.  This  approach  would  also  save 
computing  time.  Future  work  with  the  code  will  include  this  modification. 

The  code  is  currently  configured  to  model  the  experiments  conducted  at  the 
University  of  Michigan,  but  could  easily  be  altered  to  simulate  a  wide  variety  of  shock  and 
detonation  phenomena.  The  basic  FCT  subroutine  used  can  treat  rectangular,  cylindrical, 
or  spherical  one-dimensional  systems.  The  operator-splitting  method  used  here  for  the 
two-dimensional  simulation  could  also  easily  be  extended  to  three  dimensions.  It  is  hoped 
that  replacement  of  the  perfect-gas  equation-of-state  (EOS)  with  a  more  general  EOS  will 
allow  simulation  of  detonation  in  condensed  explosives.  This  approach  has  already  been 
used  by  Guirguis,  Oran  and  Kailasanath,  who  used  the  HOM  EOS  [22j  and  an  FCT  algorithm 
to  study  the  cellular  structure  of  liquid  nitromethane  [12,  131.  FCT  codes  are  also 
currently  being  developed  at  MRL  to  model  shock  propagation  in  metals,  and  preliminary 
results  are  encouraging  [231. 

All  the  calculations  reported  here  were  performed  on  the  NRL  Cray  X-MP/12 
computer.  The  code  requires  approximately  900,000  words  of  memory,  and  typical  run 
times  for  3,000  cycles  were  around  30  minutes. 
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SYMBOLS 

<t>  Equivalence  ratio,  i.e.  fuel  to  oxygen  ratio  for  a  given  system,  divided  by  the 

same  ratio  for  a  stoichiometric  system. 

p  total  mass  density 

p.  mass  density  for  species  i 

V  velocity  vector 

u,  v  x  and  y  components  of  velocity  vector 
P  pressure 

E  total  energy  per  unit  volume 

e  specific  internal  energy 

<j.  reaction  variable  for  species  i.  u  =  1  represents  pure  reactants,  u  =  0 

1  represents  pure  products. 

Ari  frequency  factor  for  species  i  in  reaction  rate  expression 

Eri  activation  energy  for  species  i  in  reaction  rate  expression 

r .  induction  time  for  species  i 

.  frequency  factor  for  species  i  in  induction  time  expression 

E^ .  activation  energy  for  species  i  in  induction  time  expression 

fj  fraction  of  induction  time  elapsed  for  species  i 

p.  fraction  of  species  i  within  a  cell 

Cvj  specific  heat  at  constant  volume  for  species  i 
mwj  molecular  weight  of  species  i 

AE.  energy  release  for  species  i 

Cj  sound  speed  of  species  i 

6t  hydrodynamic  time  step 

At  time  remaining  before  end  of  induction  period 

R*  universal  gas  constant 


Table  1 

Parameter  Values  tar  Hg.OjiAr  in  ratio  2:1:7 


aE 

- 

0.180  kcal/gm 

y 

= 

1.5555 

mw 

= 

31.6 

Ar 

= 

24.0  x  108  sec-1 

Er 

= 

3.0  x  104  kcal 

A 

r 

= 

5.6840  x  10"8  sec 

E 

r 

- 

1.6031  x  104  kcal 

FIGURE  1  Schematic  of  the  Computational  Domain 


FIGURE  2 


Pressure  Contour  Plots  :  Non-reactive  code. 
Initial  disturbance  on  left  hand  boundary,  dividing 
barrier  not  yet  implemented. 


1000  cycles 


2200  cycles 


1400  cycles 


2600  cycles 


FIGURE  3 


Pressure  Contour  Plots  :  Non-reactive  code. 

Initial  disturbance  on  left  hand  boundary  of  upper  chamber 


2500  cycles 


FIGURE  4 

Pressure  Contours  :  Single  Species,  Reactive 
Cycle  Numbers  700  to  2800 


FIGURE  4  continued 


FIGURE  5 

Pressure  Contours  :  Two  species,  one  reactive,  one  inert 
Cycle  Numbers  700  to  3000 


FIGURE  5  continued 
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FIGURE  5  continued 


FIGURE  5  continued 


FIGURE  6 

Pressure  Contours  :  Two  species,  with  SP1  =  SP2 
Cycle  Numbers  800  to  2000 


FIGURE  7  continued 


FIGURE  7  continued 
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The  development  of  a  two-dimensional  computer  code  to  simulate  detonation 
transfer  between  explosive  layers  is  described.  The  code  is  based  on  previous  models 
developed  at  the  Laboratory  for  Computational  Physics  and  Fluid  Dynamics  at  Naval 
Research  Laboratories,  'Washington .TJCHo  study  the  structure  of  layered  detonations  and 
the  details  of  detonation  transmission  from  one  medium  to  another.  The  code  has  been 
configured  to  simulate  experiments  conducted  at  the  University  of  Michigan  on  the  lateral 
transfer  of  detonation  and  shock  phenomena  between  different  gaseous  layers.  Preliminary 
calculations  with  the  code  show  that  the  computations  produce  many  of  the  structures  seen 
in  the  Michigan  experiments  and  also  provide  detailed  descriptions  of  the  detonation 
transmission  and  evolving  structure.  , t 
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